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Abstract 

The thermal properties of bulk copper are investigated by performing ab ini- 
tio DFT and DFPT calculations and using the quasiharmonic approximation 
for the free energy. Using both the LDA and the GGA for the exchange- 
correlation potential, we compute the temperature dependence of the lattice 
constant, coefficient of thermal expansion, bulk modulus, pressure deriva- 
tive of the bulk modulus, phonon frequencies, Griineisen parameters, and 
the electronic and phonon contributions to the specific heats at constant vol- 
ume and constant pressure. We obtain answers in closer agreement with 
experiment than those obtained from more approximate earlier treatments. 
The LDA/GGA errors in computing anharmonic quantities are significantly 
smaller than those in harmonic quantities. We argue that this should be a 
general feature, and also argue that LDA/GGA errors should increase with 
temperature. 

PACS #: 63.20.Ry, 65.40.+g, 65.70.+y, 71.15.Mb 
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I. INTRODUCTION 



In any study of the properties of metals, it is obviously crucial to include the effects of 
temperature. Thermal expansion results from the anharmonicity of the interatomic poten- 
tials, and this change in the lattice constant upon heating a metal is accompanied by changes 
in the elastic and vibrational properties. Experimental measurements of the temperature 
dependence of the lattice constant, elastic moduli, phonon frequencies, Griineisen parame- 
ters, etc. of most elemental metals have been available for a few decades now. However, it 
has become possible to calculate these thermal properties from first principles only in the 
last few years. 

There are two main issues to be resolved when trying to compute the thermal properties 
of metals: one is how to describe the interatomic interactions accurately, and the other is 
how to incorporate the effects of temperature into this description. 

To date, most computations of the thermal properties of metals have made use of 
parametrized interatomic potentials. This necessarily introduces errors, even when the 
potentials are semi-empirical and include both theoretical and experimental values in the 
fitting database. Self-consistent density functional theory (DFT) calculations provide the 
most accurate way of computing interatomic interactions from first principles. Using the 
DFT prescriptions to obtain the energies as a function of nuclear co-ordinates avoids the 
errors introduced by assuming parametrized forms of interatomic potentials. 

As for computing the effects of temperature, one possible approach is to perform molec- 
ular dynamics simulations at finite temperatures. This approach has, for example, been 
combined with empirical and semi-empirical potentials to calculate the thermal properties 
of metals. In principle, this approach can be extended by performing ab initio molecu- 
lar dynamics calculations at finite temperatures for large unit cells containing many metal 
atoms. However, the amount of computational effort required in order to obtain reliable 
thermodynamic averages makes this difficult, especially for the noble metals and transition 
metals, which contain tightly-bound valence electrons. Also, since such simulations treat 
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the ionic degrees of freedom classically, the results are not valid at very low temperatures, 
when zero-point effects are important. 

An alternative approach is to compute the vibrational free energy using the quasihar- 
monic approximation, in which anharmonic effects are included via the volume-dependence 
of phonon frequencies, which can be determined by performing ab initio calculations. Here 
too, in order to perform reliable averages, it is necessary to compute the frequencies for many 
wave- vectors in the Brillouin zone (BZ), which is computationally expensive, especially if 
the phonon frequencies are calculated using the "frozen-phonon" method. The development 
of density-functional perturbation theory (DFPT) jl] has considerably reduced the com- 
putational cost of obtaining phonon frequencies throughout the BZ, since unlike the frozen 
phonon method, this technique does not require using large supercells to access wave- vectors 
away from the zone center. 

Thus, combining ab initio DFPT calculations with a quasiharmonic treatment of the 
anharmonicity of vibrations currently offers us the most reliable yet practicable approach 
towards calculating averaged thermal properties, at least up to temperatures not too close 
to the melting point. In recent years, this combined approach has been shown to be quite 
successful in predicting the bulk thermal properties of the simple metals Al, Li and Na [Q, 
and the noble metal Ag ||. 

However, there remains one important issue that has to be decided when performing 
ab initio calculations: how to describe the exchange and correlation effects in the electron- 
electron interactions. The exact form of the exchange- correlation functional is not known, 
and one has to use various approximate schemes; the most widely used ones being the local 
density approximation (LDA) and various versions of generalized gradient approximations 
(GGAs). The GGAs are intended to be an improvement on the conventional LDA, and do 
indeed perform better in certain situations, such as transition states in chemical reactions, 
or systems containing "weak" bonds. Unfortunately, however, the GGAs do not always give 
answers that are in better agreement with experiment. 

Improving the treatment of exchange and correlation effects is the holy grail in the 
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field of electronic structure calculations, and as an aid towards achieving this goal, it is 
desirable to have a clear picture of the comparative merits of the LDA and GGA in various 
situations. It has been known for a long time now that the LDA tends to "overbind" , giving 
lattice constants that are too small, and bulk moduli, phonon frequencies and cohesive 
energies that are too large. The GGAs seem to overcorrect these errors, giving lattice 
constants that are too large. A recent study Q showed that this overcorrection is manifested 
also in the harmonic properties: the GGA gives bulk moduli and phonon frequencies that 
are systematically lower than the experimental ones. We are not aware of any detailed 
studies comparing the performance of the LDA and GGA in describing anharmonic effects, 
which manifest themselves in the temperature-dependence of the lattice constant, elastic 
and vibrational properties, and specific heat capacities, and in the values of anharmonic 
quantities such as the Griineisen parameters. 

To this end, in this paper, we have performed ab initio calculations to study the ther- 
mal properties of bulk copper, using both the LDA and GGA. We have computed the 
temperature-dependence of the lattice constant, the coefficient of thermal expansion, the 
isothermal bulk modulus, the phonon frequencies, the individual and overall Griineisen pa- 
rameters, and the specific heat capacities at constant volume and constant pressure. 

II. AB INITIO CALCULATIONS 

The ab initio calculations were performed using the PWSCF and PHONON codes ||. 
Total energies were computed using DFT, and phonon frequencies using DPFT. The interac- 
tion between the ions and valence electrons was described using an ultrasoft pseudopotential 
0. A plane- wave basis set with a cut-off of 30 Ry was used; a cut-off of 300 Ry was used 
in the expansion of the augmentation charges necessitated by the use of the ultrasoft (non 
norm-conserving) pseudodopotential. Brillouin-zone integrations were performed using 60 k 
points in the irreducible part of the BZ. Phonon dynamical matrices were computed ab initio 
for a 4 x 4 x 4 q point mesh; Fourier interpolation was then used to obtain the dynamical 
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matrices on a 24 x 24 x 24 q point mesh. This latter set was used to evaluate all quantities 
that involve an integration over phonon wave-vectors q. 

In order to deal with the possible convergence problems for metals, a smearing technique 
was employed using the Methfessel-Paxton (MP) scheme 0, with the smearing parameter 
a set equal to 0.05 Ry. However, when evaluating the electronic contribution to the specific 
heat capacity (as explained below) we instead used a Fermi-Dirac (FD) smearing, with the 
electronic levels occupied according to the FD distribution appropriate to the temperature 
of interest T. (Incidentally, with this latter scheme, we did not face convergence problems 
even at low values of T.) 

When using the LDA, we used the parametrization by Perdew and Zunger of the results 
of Ceperley and Alder ||. For the GGA, we used the Perdew-Burke-Ernzerhof form ||; this 
choice was made in part because it is easier to implement in the DPFT calculations, and 
because it gives a good description of the linear response of the uniform electron gas. 

To summarize, the following were obtained from DFT and DPFT calculations: (i) Total 
energies at a range of lattice constants, using (a) MP smearing (b) FD smearing, for a range 
of temperatures between 1 K and 1400 K. (ii) For each lattice constant, the dynamical 
matrices (and thus phonon frequencies) for the 4 x 4 x 4 set of q points, using MP smearing; 
Fourier interpolation was then used to obtain the dynamical matrices on the 24 x 24 x 24 
set of q points. (It was verified that replacing the MP smearing by the FD smearing did not 
make an appreciable difference to the phonon frequencies, i.e., the latter are not sensitive 
to the electronic temperature.) All of the above quantities were computed using both the 
LDA and the GGA. This set of results was then used to calculate the thermal behavior, as 
described below. 

III. RESULTS AND ANALYSIS 

The static results for lattice constant do, the bulk modulus Bq and the pressure-derivative 
of the bulk modulus B' are obtained by fitting the results for the static total energies (using 
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MP smearing) versus lattice constant to the fourth order Birch-Murnaghan equation of 
state 0. Using the LDA, we obtain a = 6.71 bohr, B = 1.72 MBar, and B' = 5.0. 
The corresponding results with the GGA are a® = 6.94 bohr, Bq = 1.28 MBar, and B' = 
5.11. As expected, the experimental values for the lattice constant (ao = 6.82 bohr) (Tl| 
and bulk modulus (1.37 Mbar) [|I^] lie sandwiched between the LDA and GGA values; it 
should however be noted that the experimental values are at room temperature, and the 
calculated values listed above do not yet include the effects of temperature. For B', there 
does not seem to be a consensus on the experimental value, with several values reported in 
the literature. Listed in chronological order, these are: 3.91 fll3 |, 5.3 [fL4|, 4.8 [15[], 4.1 ||16|| , 
5.59 |7| and 5.44 @. 

To study the effects of changing temperature, one has to look at the free energy, incor- 
porating the effects of thermal vibrations (phonons). The free energy at temperature T and 
lattice constant a is given, within the quasiharmonic approximation, by: 

qA 

Here, the first term on the right-hand-side is the static energy E stat (a), and the second term 
is the vibrational free energy. The sum is over all three phonon branches A and over all 
wave- vectors q in the BZ (we will use the 24 x 24 x 24 q mesh in evaluating this), h is 
Planck's constant, ks is Boltzmann's constant, and u}q\(a) is the frequency of the phonon 
with wave- vector q and polarization A, evaluated at lattice constant a. 

The lattice constant at temperature T, ao(T) is obtained by minimizing F(a,T) with 
respect to a. The linear expansion e(T) is then given by 

<t) = a ° (r) -;° (r ' ) , ( 2) 

where T c is the reference temperature of 298.15 K. 

Fig. 1 shows the results for e(T), (expressed as a percentage) using both the LDA and 
the GGA, compared to the experimental value [|BJ. It is seen that the agreement with 
experiment is quite good, though the LDA slightly underestimates the expansion, and the 
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F{a,T) = E stat (a) + k B Tj2ln^2smh(^- 



GGA slightly overestimates it. This becomes more obvious upon differentiating the results 
for a Q (T) to obtain the coefficient of linear expansion: 

(Note that this definition of ot(T) is the one used for experimental data. When using a(T) 
in thermodynamic relations, ao(T c ) should be replaced by clq(T) in the right-hand-side of 



the above equation.) Fig. 2 compares the calculated and experimental values |19| for a(T) 
up to a temperature of 1400 K (the experimental value for the bulk melting temperature 
is 1357 K). Once again, it is clear that the experimental values lie sandwiched between the 
LDA and GGA values, though they lie somewhat closer to the LDA values, especially at 
high temperature. However, it should be pointed out that the calculated values may be 
inaccurate at very high temperatures for two reasons: (i) The use of the quasiharmonic 
approximation may not be justified at temperatures just below the melting point, as this is 
expected to be a region of high anharmonicity. (ii) A part of the experimentally measured 
thermal expansion at high temperatures results from the formation of vacancies; this effect 
is not included in our calculations, where we assume that the crystal remains defect-free at 
all temperatures. 

By fitting the results for the free energy from Equation 1 to the fourth-order Birch- 



Murnaghan equation of state flQ |, we also obtain the variation in temperature of the bulk 
modulus B and the pressure derivative of the bulk modulus B'. We find that the quality of 
the fit is noticeably better with the fourth-order equation of state than with the Murnaghan 



equation of state |2(| or with the third-order Birch-Murnaghan equation, especially at higher 



temperatures. The results for Bq(T) are plotted in Fig. 3, from which it can be seen that 
though, at all temperatures, the absolute value of B (T) is overestimated by the LDA and 
underestimated by the GGA, the rate of change of B with temperature is approximately 
the same for both, and moreover, this rate agrees well with that measured experimentally 
|nj . Fig. 4 shows the results for B'(T); it can be seen that B' depends noticeably on the 
temperature. (Incidentally this temperature dependence is considerably underestimated if 



one uses the Murnaghan equation or the third-order Birch- Murnaghan equation.) At 300 K, 
the LDA and GGA values for B' are 5.21 and 5.40, respectively, compared to the static values 
of 5.00 and 5.11. These values agree well with some of the room-temperature experimental 



values cited above p3|-fi8[, though there is a considerable scatter in the experimentally 
reported values. 

Since we know how the phonon frequencies vary with ao, and how ao varies with T, it is 
now a simple matter to get the phonon frequencies at any desired temperature. In Fig. 5, 
the calculated and measured |22j phonon frequencies, at a temperature T=80 K, are plotted 
along several high-symmetry directions in the Brillouin zone. At this temperature, the LDA 
gives a = 6.73 bohr, and the GGA gives a = 6.96 bohr (since the temperature is relatively 
low, there is not an appreciable change from the static values). Yet again, it can be seen that 
the experimental values lie in between the LDA and GGA values. The over/under-estimation 
of the frequencies by the LDA/GGA can be traced back to the under/over-estimation of 
the lattice constant. In fact, if the phonon frequencies are computed at the experimental 
lattice constant, the situation is reversed, and the GGA frequencies are higher and the LDA 
frequencies are lower than experiment, though the latter are closer to the experimental 
values than the former. 

We can also compute the temperature-dependence of the specific heat capacities at con- 
stant volume and constant pressure, as described below. The specific heat at constant 
volume has two contributions, one from the phonons, and the other from the electrons. The 
former is given by 

' au; qA (a (T))\ 



2 i 



(4) 



qA qA 

yell 



2k B T J 

The electronic contribution to the specific heat, Cy(T) is obtained from the self-consistent 
DFT calculations using FD smearing corresponding to a temperature T, by computing the 
derivative with respect to the smearing temperature T of the electronic entropy, evaluated 
at the corresponding lattice constant cto(T). The total specific heat at constant volume is 
then C^{T) = Cf(T) + CfJ(T). 
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C p , the specific heat at constant pressure, can then be computed by using the relation: 
C P (T) = C$\T) + ^a 2 (T)B (T)a (T)T. (5) 

Figs. 6 (a) and (b) show the results thus obtained for C£ h (T), C$(T), C^(T) and C P (T), 
computed using the LDA and GGA respectively. As expected, the electronic contribution to 
the specific heat, Cy(T), is much smaller than the phonon contribution Cy h (T), though not 
negligible. The experimental values for C P (T) Jl9[ are also plotted. It is seen that for both 
the LDA and GGA, the agreement with experiment is excellent up to about 600 K. Above 
this temperature, the agreement remains very good for the LDA, but is poorer for the GGA. 
Note that at these high temperatures, Cy- h (T) has reached its saturation value of 3k b per 
atom, and the LDA and GGA values for Cy h (T) are therefore identical. For the difference 
between C p and Cy, the error due to the under /over-estimation of a by the LDA/GGA is 
to some extent cancelled out by the over/under-estimation of Bq. 

The anharmonicity of the vibrations can be examined by computing the mode Griineisen 
parameters, defined by 

where V = a 3 /4 is the volume of the unit cell. Fig. 7 shows the results for the Griineisen 
parameters for the same high-symmetry modes for which the frequencies were plotted in 
Fig. 5. They have been evaluated at the static lattice constants. Though the LDA and 
GGA static lattice constants are different, it can be seen that the discrepancy in the corre- 
sponding Griineisen parameters is small; considerably smaller than the differences in phonon 
frequencies. For example, at the X point (zone edge along [100]), the discrepancy between 
the LDA and GGA results for the phonon frequencies is 10.7% and 12.5% for the transverse 
and longitudinal branches respectively, whereas the corresponding Griinesien parameters 
differ by only 0.4% and 2.1% respectively. It is interesting to compare Fig. 7 with Fig. 3 
of Ref. H, which shows the corresponding result for Ag. Though the phonon dispersion 
curves of Ag and Cu are very similar in shape and structure, our curves for the Griineisen 



parameters of Cu look quite different, in some areas of the BZ, from those reported earlier 
for Ag. 

We also compute the overall Griineisen parameter 7, which is obtained by averaging over 
the individual Griineisen parameters 7 q A of all the modes, using the equation: 



where the contribution from each mode (qA) is weighted by CV(qA), its contribution to the 
specific heat, as defined in Equation 4. This quantity is of interest because it appears in 
some useful thermodynamic relations (as discussed below) and experimental papers often 
report this overall value. The temperature dependence of 7 comes from the temperature 
dependence of both the individual Griineisen parameters 7 q A (which depend on the lattice 
constant and hence on T) and that of the specific heat. 

Fig. 8 shows the results for the results for the variation of the overall Griineisen parameter 
7 as a function of temperature. It is seen that the percentage difference between the LDA 
and GGA results is small, ranging from 2.5% at 100 K to 5.3% at 1300 K. The discrepancy 
with experiment is also quite small, with the LDA and GGA errors being 2.8% and 4.8% 
respectively at room temperature. In comparison, the LDA and GGA results for the bulk 
modulus (plotted in Fig. 3) differ from each other by 30 to 50% over the same temperature 
range, with the LDA and GGA errors (with respect to experiment) being 18.8% and 13.7% 
respectively at room temperature. 



We compare our results to those of three previous calculations. In the first, interatomic 
potentials are described by a pair potential fit to experimental data, and thermal effects are 
treated by formulae that are valid in the high-temperature limit. In the second, the static 
energies at zero temperature are computed ab initio; however, thermal effects are computed 
using a Debye model and various approximate relations; i.e., the phonon frequencies are not 
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E q ACV(qA) 



IV. COMPARISON WITH EARLIER CALCULATIONS 



10 



calculated ab initio. In the third calculation, the interatomic potentials are described by an 
empirical form that includes some many body effects, and thermal effects are treated using 
finite-temperature molecular dynamics simulations. 

In their study of the thermodynamic properties of face-centered-cubic (FCC) metals, 
MacDonald and MacDonald [^J have described the interatomic interactions using a modified 
Morse potential fit to experimental data such as the Debye temperature G^, the sublimation 
energy, and the thermal expansion in the neighborhood of Bd (342 K). The electronic 
contribution to Cy is estimated using free-electron theory. Though the thermal expansion is 
fit to agree with experiment at low temperatures, a is underestimated by about 20% at 1200 
K (in comparison, our results for a at 1200 K, with no fit to data on thermal expansion, 
differ from experiment by -12% when using the LDA and +29% when using the GGA.) The 
absolute value of B Q agrees well with experiment (which is to be expected, given the fitting to 
Qd)', however 3Bq{T) / dT is underestimated. The magnitude of the electronic contribution 
to Cy, estimated from free-electron theory, is similar to what we obtain from our more exact 
approach; and the agreement between the calculated and experimental values for C p is fairly 
good, similar to that obtained by us. The value of 7 changes from 1.947 at 100 K to 2.127 
at 1000 K, which is comparable to our results. However, it should of course be kept in mind 
that in our calculations we do not fit to any empirical data at all. 

Moruzzi, Janak and Schwarz |24| have computed the bulk binding curve (i.e., E stllt (a)) 
by performing ab initio augmented-spherical-wave method calculations. They use this to 
obtain B (and thus an approximate Op) and 7 (independent of T). The free energy is 
then evaluated in the Debye model, with the volume-dependence of the frequencies being 
determined by 7. They evaluate a only up to T = 300 K, getting a value of a that is too 
low at 300 K by 20%. Our calculation improves upon this one in that we do not use a Debye 
model, and do not assume that all modes have the same degree of anharmonicity, since we 
calculate individually and exactly the values of ujq\(a). This is probably why our calculated 
values for a are closer to experiment: at 300 K, our errors in the calculated a are -14% 
(LDA) and +14% (GGA). 
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Cagin et al. |2l| have used the empirical Sutton-Chen potential to describe the inter- 
atomic interactions. The potential parameters are fit to the cohesive energy, bulk modulus, 
etc., at K. Temperature effects are determined by performing molecular dynamics simu- 
lations. As a comparison, let us consider the values for the thermal expansion e at 1000 K: 
they obtain a result of 2.42%, compared to our values of 1.22% (LDA) and 1.64% (GGA), 
and the experimental value of 1.37%. Our values are clearly closer to experiment; in this 
case, their larger errors presumably arise from deficiencies in their interatomic potential. 

To summarize: though our results do not agree exactly with experiment, we still do a 
better job than earlier calculations. This is because we have eliminated the errors due to 
the utilization of parametrized interatomic potentials and an approximate treatment of the 
lattice vibrations. The (smaller) errors that remain in our calculations are due to the choice 
of exchange-correlation potential (and, possibly, the use of the quasiharmonic approximation, 
though we believe these errors to be small). 

V. DISCUSSION OF RESULTS 

From the results presented in Section III, it is clear that the approximations used for the 
exchange-correlation potentials introduce much smaller errors in anharmonic quantities such 
as 7, B' and dB (T)/dT, than in harmonic properties like the bulk modulus and phonon 
frequencies. 

At first sight, the relatively large discrepancy between the LDA and GGA values for the 
coefficient of thermal expansion a may seem to contradict this statement. However, a is 
not a purely anharmonic quantity. For example, in a one-dimensional anharmonic potential 
given by V(x) = \cx 2 — Igx 3 , we have a oc g/c 2 [0, i.e., a depends on both the harmonic 



coefficient c and the anharmonic coefficient g, and an error in the former will be manifested 
as an even larger error in a. For the present case, where we have to average over a number 



of normal modes in three dimensions, it is possible to derive [11| a corresponding equation 



relating a to the averaged anharmonic quantity 7 and harmonic quantity B : 
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a(T) 



l(T)C v (T) 
3B (T) 



(8) 



From Fig. 8, it is seen that the discrepancy between the LDA and GGA values of j{T) 
is small, and that both are close to experiment. Any error in CV(T) is negligible, especially 
at temperatures above the Debye temperature, where Cy has reached its saturation value 
of ?>k B per atom. There remains the large error in B (T). At a temperature of 1000 K, for 
example, the LDA and GGA values for 7, B and a differ by 4%, 41%, and 32% respectively, 
which is consistent with our argument that the discrepancy in a arises almost entirely from 
the discrepancy in B . 

If we assume that the main source of the LDA and GGA errors in the values of physical 
quantities is the wrong value obtained for the lattice constant, then it is indeed consistent 
that the errors in anharmonic quantities should be smaller than those in harmonic quantities. 
For example, let us return to the example of the one-dimensional cubic anharmonic potential 
cited above: 



If the second derivative (which gives harmonic properties) is evaluated at a lattice constant 
a which is not the correct lattice constant a , the error in the second derivative is given by 
d 2 V/da 2 \ a —c = —g(a — a ), and the greater the error in the lattice constant, the greater 
the error in the calculated harmonic properties. However, there is no error in the third 
derivative (which gives cubic anharmonic properties), since d 3 V/da 3 \ a — —g for all a. (If 
one were to add a quartic term to Eq. 9, then the error in the second derivative would have 
terms linear and quadratic in (a — ao), whereas the error in the third derivative would be 
proportional to (a — a ), and the fourth derivative would be exact.) Thus, errors arising 
from using a wrong lattice constant are manifested to lesser and lesser degrees as one goes 
to higher order anharmonic properties. 

One can also argue that the LDA and GGA errors made in computing physical properties 
should increase with temperature: At T=0, the LDA underestimates ao an d the GGA 




(9) 
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overestimates it. Upon heating, the LDA underestimates the thermal expansion (because 
of the overestimation of B ) and the GGA overestimates it (because of the underestimation 
of Bq). Thus, as the temperature is increased, the underestimation by the LDA and the 
overestimation by the GGA of the lattice constant are both aggravated further, resulting 
in increasingly unreliable results for physical properties, though the errors are smaller for 
anharmonic properties than harmonic ones. To avoid this, we suggest that when performing 
ab initio DFT calculations at high temperatures in systems that display large errors in the 
calculated bulk modulus, it is perhaps a good idea to use the experimental values of the 
thermal expansion, regardless of whether one is using the LDA or GGA. It is also useful to 
keep in mind that the error in the static value for the bulk modulus is already a good indicator 
of the magnitude of the error that will be made in the coefficient of thermal expansion; thus, 
if in a particular case, either the LDA or the GGA gives a better value for ths static value 
of ao and Bq, it will probably also give a better description of finite-temperature properties. 



VI. SUMMARY 

To summarize: we have performed ab initio calculations to study the thermal properties 
of bulk copper, using phonon frequencies computed using DFPT, and the quasiharmonic ap- 
proximation for the vibrational free energy. We have calculated the temperature dependence 
of the lattice constant do, the thermal expansion e, the coefficient of thermal expansion a, 
the bulk modulus B , the pressure derivative of the bulk modulus B', the overall Griineisen 
parameter 7, and the lattice and electronic contributions to the specific heat capacities and 
constant volume and pressure, Cy and C p . We have also presented results for the dispersion, 
along high-symmetry directions in the BZ, of the phonon frequencies and mode Griineisen 
parameters. All of the above have been computed using both the local density approximation 
(LDA) and generalized gradient approximation (GGA). 

Neither the LDA nor the GGA is clearly to be preferred in this case, with both giving 
errors of comparable magnitude (though generally of opposite sign). At all temperatures, 
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the LDA systematically underestimates the lattice constant and the coefficient of thermal 
expansion, and the GGA overestimates these. In contrast, the LDA always overestimates the 
bulk modulus and phonon frequencies, and the GGA underestimates them. The electronic 
contribution to the specific heat is found to be considerably smaller than the phononic 
contribution, though not neglible. The results for C p agree very well with experiment, 
except for the GGA results at high temperatures. However, the discrepancy between the 
LDA and GGA results (and their discrepancy with experiment) is considerably lower for 
anharmonic quantities such as 7, B' and dB (T)/dT than for the harmonic properties. In 
any event, our results are closer to experiment than those of earlier calculations in which 
the interatomic interactions and/or thermal effects were treated approximately. 

We have argued that if the main source of errors can be attributed to the wrong value 
obtained for the lattice constant resulting from the approximate nature of the exchange- 
correlation potential, then it is indeed reasonable that the errors in anharmonic quanti- 
ties should be smaller than those in harmonic quantities. We have also argued that the 
LDA/GGA errors should increase with temperature, suggesting the need for care and cau- 
tion when performing ab initio calculations at high temperatures. 
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FIGURES 

FIG. 1. Linear thermal expansion e as a function of temperature, referred to a reference tem- 



perature T c of 298.15 K. The experimental values are from Ref. [19|. It is seen that both the LDA 



(solid line) and GGA (dashed line) results are close to the experimental values. 

FIG. 2. Coefficient of linear thermal expansion a as a function of temperature. Both the LDA 
and GGA values are reasonably close to the experimental values; however the LDA underestimates 



and the GGA overestimates the thermal expansion. Experimental values are from Ref. [19]. 



FIG. 3. Variation with temperature of the bulk modulus Bq. At all temperatures, the LDA 
(solid line) overestimates -Bo and the GGA (dashed line) underestimates it; however, dB^/dT is 
approximately the same for the LDA, the GGA, and the experimental values. The experimental 



values are from Ref. |21| . 



FIG. 4. Variation with temperature of the pressure derivative of the bulk modulus, B'. The 
solid line is the LDA result, and the dashed line is the GGA result. 

FIG. 5. Phonon dispersion along high-symmetry directions in the BZ, at 80 K. The solid and 
dashed lines are the results obtained using the LDA and GGA respectively, and the filled circles 



are the experimental values from Ref. [22|. "L" and "T" denote the longitudinal and transverse 
branches respectively. 

FIG. 6. Calculated values of C p and C v , in units of ks per atom, obtained using the (a) LDA 
and (b) GGA. The dot-dashed lines show Cy 1 , the thin dashed lines show Cy , and the thick dashed 
lines show their sum Cy*. The solid lines show the calculated values for C p , obtained from Cy* 
by using Eq. (5). The dots show the experimental results for C p , as given in Ref. ||. 

FIG. 7. Calculated dispersion curves for the individual mode Griineisen parameters 7 q ^ for 
the same high-symmetry directions in the BZ as in Fig. 5. The values have been evaluated at 
the static lattice constant; the discrepancy between the LDA (solid lines) and GGA (dashed lines) 
results is small. "L" and "T" denote the longitudinal and transverse branches respectively. 
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FIG. 8. Overall Griineisen parameter 7 as a function of temperature. The solid and dashed 
lines are the results obtained using the LDA and GGA respectively. The experimental value is 



taken from Ref. 24]. 
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